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We introduce a general technique to compute finite temperature electronic properties by a novel 
covariant formulation of the electronic partition function. By using a rigorous variational upper 
bound to the free energy we are led to the evaluation of a partition function that can be computed 
stochastically by sampling electronic wave functions and atomic positions (assumed classical). In 
order to achieve this target we show that it is extremely important to consider the non trivial 
geometry of the space defined by the wave function ansatz. The method can be extended to any 
technique capable to provide an energy value over a given wave function ansatz depending on several 
variational parameters and atomic positions. In particular we can take into account electronic 
correlation, by using the standard variational quantum Monte Carlo method, that has been so far 
limited to zero temperature ground state properties. We show that our approximation reduces 
correctly to the standard Born-Oppenheimer (BO) one at zero temperature and to the correct high 
temperature limit. At large enough temperatures this method allows to improve the BO, providing 
lower values of the electronic free energy, because within this method it is possible to take into 
account the electron entropy. We test this new method on the simple hydrogen molecule, where 
at low temperature we recover the correct BO low temperature limit. Moreover, we show that the 
dissociation of the molecule is possible at a temperature much smaller than the BO prediction. 
Several extension of the proposed technique are also discussed, as for instance the calculation of 
critical (magnetic, superconducting) temperatures, or transition rates in chemical reactions. 



I. INTRODUCTION 



The calculation of finite temperature electronic properties is one of the most important and challenging aspects of the 
numerical simulations. In the past several progress have been done by extending the DFT method to finite temperature 
1, 2] or by using quantum Monte Carlo yl (QMC) within various path integral formulations especially in the 

study of the hydrogen phase diagram [1CH17I ]. In both cases many problems remain as for instance the lack of an 
accurate local functional at finite temperature for DFT methods prevents so far practical applications, and, within 
QMC techniques, the difficulty to deal with the fermion sign problem[3], restricts the spectrum of applicability to 
very few cases and very limited temperature ranges. On the other hand it is clear that, in many physical phenomena, 
such as the occurrence of magnetic or insulating phases below a critical temperature, the electronic entropy cannot 
be neglected, even when the small ratio A e i between the electronic mass and the atomic one, allows the decoupling of 
the electronic degrees of freedom from the atomic ones, within an acceptable approximation. In this paper we aim to 
extend the validity of the Born-Oppenheimer approximation in the following sense. By using the smallness of X e i we 
are generally lead to compute an electronic partition function Z\R] at fixed nuclei position: 

Z = J dR Z\R] (1) 
Z[R] = Trexp(-ff R /T) (2) 

where T is the temperature (here and henceforth the Boltzman constant is assumed to be one and we neglet for 
simplicity the overall constant coming from integration of the atomic momenta), i?R is the standard electronic 
Hamiltonian, that includes also the classical ionic contribution, and that depends only parametrically upon the 
atomic positions R. Eq.([T]) is the first step of the Born-Oppenheimer approximation that- we remark- is generally 
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valid for A e i small, namely when the temperature T is large enough that quantum effects on protons can be neglected. 
The second approximation, usually adopted within the BO approximation, is to assume that the electronic degrees 
of freedom have a gap much larger than the temperature T so that Z[R] can be approximated by exp(— Eq(R)/T) 
where Eq(R) is the ground state energy of the hamiltonian i?R. In the following derivation we want to avoid the latter 
approximation, because, as emphasized before, in several cases it may fail even when we are in the limit of small X e i. 
For instance the occurrence of a broken symmetry phase often implies gapless electronic excitations in Hr, and the 
approximation Z[R] — exp(— Eo(R)/T) cannot be safely assumed. Other examples are conical intersections [l9l - [2l| . 
when for some particular ionic positions i?R becomes gapless and nearby the proximity between different (namely 
corresponding to low-lying excited states) BO energy surfaces is possible. In this conditions a pure electronic ground 
state technique fails as the tunneling between different BO energy surfaces cannot be taken into account consistently. 
As the last very important example we mention the calculation of transition rates in chemical reactions, that cannot 
be accurately computed within a pure BO approximation [22l - (26| . 

The main task of this paper is to device a method, able to quantify finite temperature properties of realistic systems, 
within a rigorous variational upper bound of the total free energy F — — TlnZ, in the limit of small X e i. The method 
we propose is supposed to be simple enough to avoid most of the known drawbacks, as does not rely on the knowledge 
of any particular functional, or, within our variational approximations, can be employed by quantum Monte Carlo, 
without facing the so called "fermion sign problem" . 

The paper is organized as follows. The derivation of the approximate expression of the electronic partition function 
introduced and used in this work is given in section UH and some important but more detailed aspects are reported 
in appendixes E] [B] [C] and [D] This derivation is not specific for a QMC framework, indeed App. [D] is specifically 
oriented to an implementation of the method into a Hartree-Fock or DFT framework. Next we show how to sample 
the introduced partition function using a Langevin dynamics for the wave function parameters, sections IIII[ and the 
ion coordinates, section IPV1 In section [V] we finally show some results we have obtained using this approach for the 
hydrogen molecule. 



II. FINITE TEMPERATURE ELECTRONIC PARTION FUNCTION 



We consider the problem to estimate the finite temperature partition function of an electronic system with N 
electrons and M atoms, where we assume in the following that, as discussed in the introduction, the ions are classical 
particles, whose coordinates R appear just as simple parameters in the electronic hamiltonian H^i and are confined 
in a finite volume V. Therefore, once the ion positions are fixed, we need to evaluate the electronic partition function: 



Z[R] = Trexp(- ( 3if R ) 



(3) 



where f3 — 1/T. Our derivation applies for an Hamiltonian with a bounded spectrum defined in a finite Hilbert 
space with dimension D. For instance in electronic structure calculation one can consider a finite dimensional basis of 
localized orbitals around each atom. In order to simplify the notations we can consider standard creation operators 
with canonical commutation rules, spanning the finite single electron basis, as for a standard lattice Hamiltonian, 
namely cj for i = 1, • • ■ ,L, where for shorthand notations i labels also the spin, namely i < L/2 (i > L/2) refers to 



spin up(down)-states. We consider the generic wavefunction 



J x \SD), where: 



J = exp(l/2^ 



\SD) 



N L 

i=i j=i 



ct 



10) 



(4) 
(5) 



and rii = cjci, for a system of N electrons. In the continuous limit this wave function is the standard Jastrow-Slater 
one used in quantum Monte Carlo in order to describe electron correlation [27] . Estensions of this wave function are 
possible using AGP 28], Pfaffian 29], backflow[30], and the following considerations apply also for these more recent 
ansatz, because they all contain the Slater determinant \SD) in a particular limit. 

In all cases, the real variational parameters that define the above wave function, namely Vij and ip^ are compactly 
denoted by a = {oLi\i = \^,^ p and since all physical quantities do not depend on the norm of the wave function, we 
consider the a— manifold of states: 



(0) 
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The metric in this manifold becomes non trivial as, by a straightforward calculation, the distance between two 
states | a) and \a + da) is given by: 

ds 2 = \\\a + da) - \a)\\ 2 = -/n'r/n'.S',., (7) 

where summation over repeated indices is assumed, and S is a p x p matrix defining the metric tensor of this rather 
non trivial space, described by p independent variational parameters (e.g. a subset of Vij and ipij). The matrix S 
can be explicitly evaluated and depends only on average first derivatives of the wave function with respect to the 
parameters a' s: 



It defines a metric as it is strictly positive definite if all p variational parameters are independent and therefore 
its determinant \S\ is non vanishing. This matrix turns out to be exactly the one used in several optimization 
techniques |3ll. [32j. and can be computed also for correlated systems by sampling the correlations of the quantities 
Oj(x) = over the configuration space {x} where electrons have a definite spin and positions, namely: 

Sij = (O i O j )-(O i )(O j ) (9) 

where the symbol (. . .) denotes average over a distribution TL(x) cx (x\ip) 2 , that can be sampled by standard variational 
Monte Carlo. 

In Eq.@, we use a simple relation for recasting the trace in a finite dimensional Hilbert space as an integral of 

D 

normalized wavefunctions \c) — x i\i)i namely: 

i=l 

D J dx D 5{\\x\\ - l)(c\exp(~(3H)\c) = S D Trexp(-/3H) (10) 

where S D = 2tt d ^ 2 /T(D/2) is the area of the D— dimensional unit sphere. We note that this simple relation can be 
used to establish within a rigorous mathematical framework the finite temperature Lanczos method used in Refl33l. 
In this technique finite temperature estimates of the partition function are obtained with a finite set of randomly 
generated states |c), once it is assumed that (c\ exp(— /3H)\c), can be computed with high accuracy with the Lanczos 
method. Indeed this is nothing but evaluating statistically the integral in the LHS of Eq. fTU)) . and one does not 
need any further assumption to validate the method, apart from the fact that error bars have to be computed with 
standard statistical techniques. 

The simple relation (jlOp can be also extended in the space a with non trivial metric, by using the invariant measure 
da p y/\S\, corresponding to the metric tensor S: 



J d g*VW\(«\exp(-PHn)\*) =Trexp( _^ R ) 
Zs 

This relation is proven in App[3] provided the dimension of the space is large enough, namely contains at least the full 
space of Slater determinant wave functions, where the overall constant has been obtained by using that Z{R] = D for 



P = 0, as the metric normalization Z$ is defined as Zs — da -q — ■ We emphasize here that the relation (TTT1) is exact 
even when the dimension of the space p is much smaller than the dimension of the Hilbert space. For instance for 
real Slater determinants the number p < NL as they are defined by N orbitals each depending on L coefficients (see 
EqHJ, whereas the Hilbert space dimension D grows exponentially with L and N (See App[D]for the parametrization 
of an arbitrary real Slater Determinant)) 

In practice the number p of variational parameters defining the wave function ansatz can be much smaller than 
that necessary to span all possible Slater determinants. In the case p <C NL we expect that the equation (fTTj) is 
still valid but the trace in the RHS is limited to the largest subspace with dimension D s spanned by the variational 
ansatz. Moreover a weak dependence on R in Zs is also expected when a basis dependent on the atomic positions is 
used (it is not the case for a plane wave basis for instance). The calculation can be meaningful also in this case after 
a careful study of the dependence of the results upon the dimension of the basis chosen, as it is common practice 
in quantum chemistry calculations. In fact, in the limiting case when the one particle basis set used to define the 
orbitals in the Slater Determinant becomes complete the metric normalization Zs is independent of R, because any 
change of basis is equivalent in this limit to a mapping a — > a' . Thus Zs, being explicitly covariant, is independent 
of R and can be considered as an irrelevant constant. Therefore, within the completeness assumption, following the 
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simple derivation of App|B] we can easily bound the exact electronic partition function Z[R], because, due to the 
convexity of the exponential function, the expectation value of an exponential operator over a normalized state \a) 
satisfies: 

(a\exp(-pH R )\a) > exp(-/3(a\H n \a}). 
This immediately provides a rigorous lower bound Zq for the partition function Z: 



7>7 JrfR/rfa p ^exp(-/3( a |H R |q)) 

Z > Zq = (12) 

and a corresponding upper bond Fq for the free energy F = —TlaZ : 

F <F Q = -TlnZ Q (13) 

In this way it is evident that Fq represents an improvement to the standard Born-Oppenheimer (BO) approximation. 
In fact in this approximation only one state is assumed to contribute to the integral in Eq. (|12[) . namely the lowest 
energy state of i?R within the ansatz given by \a): 

E B o[R}=mm{(a\H R \a)} (14) 

a 

Indeed it is clear that F = minR {£?bo[R]} only at T = 0, and represents a very bad approximation to F as long as 
the temperature is raised, whereas the approximate partition function Fq approaches the correct large temperature 
limit —T\n(DV M ) of the exact partition function, while remaining a rigorous upper bound for any T. 

In App[C] we see in detail a comparison between the approximated partition function Zq here introduced, and the 
exact and BO ones, showing that our approximation turns out to be better than the BO one, above a temperature 
T* , that remains meaningful in the thermodynamic limit. 



III. MONTE CARLO SAMPLING OF THE PARTITION FUNCTION Z Q 

In principle the partition function Zq can be sampled by almost standard Monte Carlo methods, whenever the 
metric S and the expectation value of the energy H over the ansatz \a) are known, for instance within the Hartree- 
Fock theory, namely when \a) represents just a simple Slater determinant. It is also possible to replace in Zq the 
expectation value of the energy with any DFT functional depending on |a), through the corresponding density or 
gradient, the condition of functional minimum being recovered correctly at T = 0. For a discussion about the space 
of parameters for a Slater determinant wave function, and the introduction of an invariant measure in this space see 

Appini 

However in the truly correlated case, namely when the ansatz \a) differs from a Slater determinant, there are extra 
complications because both the matrix S and (a\H-£i\a) are known only within statistical accuracy. In this case a 
possible way to sample the partition function Zq and corresponding thermodynamic quantities is to use the penalty 
method[33|. introduced some years ago, by using a cost function 

V P (a, R) = (a\H n \a) - ^ In \S\ (15) 

that can be computed statistically with corresponding error bars. 

In the following we have chosen a different route, by employing a finite temperature molecular dynamics rather than 
Monte Carlo sampling, because recent quantum Monte Carlo packages provide efficient estimates of energy derivatives 
and ionic forces [35j. 

Our goal is to sample points in the electronic parameter space a distributed according to the probability distribution 
defined in Eq. (|12[) . by using first order derivatives of the cost function. In the standard Cartesian metric it is common 
practice to use a Langevin dynamics for the variables {a} and {R}, by means of the standard first order equation of 
motions (unit mass is assumed for simplicity) f36jj : 

£ = -d s V + ij (16) 

where x is a covariant vector in a finite dimensional euclidean space, whereas dgV(x) is the derivative (force) of a 
potential V. By means of this equation it is well known that it is possible to sample the equilibrium distribution 
W eq (x) = exp(— f3V(x)) provided we satisfy the fluctuation dissipation theorem given by: 



(17) 
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Now we suppose to change the reference coordinate system by means of a generic transformation of variables x — > a 
(an N— dimensional non linear mapping as in general relativity). Be the Jacobian of the transform given by the 
matrix L: 

Li.j = d Xj oti(x) (18) 
The Langevin equation in this new reference can be easily obtained: 



, dV 

d = -S- 1 ^ + Lff (19) 
aa 

where S -1 = LL\ and the equation (|17p that defines the fluctuation dissipation theorem remains unchanged. 

The Eq. (|19l) is covariant if we just replace the matrix S with the matrix defining the metric in a generic curved 
space: 

ds 2 — Si^dctidctj (20) 

where, as usual, in this formalism repeated indices are assumed summed. Indeed after the given transformation the 
above metric tensor transforms as: 

S^iL^-'SL- 1 (21) 



that, as it should, leaves unchanged the covariant first order Langevin equation (|19p . 

Thus, from the above equation, we obtain the desired result with the matrix L given by any solution of the matrix 
equation: 

S~ 1 =LLl 

Unfortunately Eq. (|19[) looks a bit complicated when it is discretized in times t n — An, because the integral of 
the random noise depends explicitly on the curvature of the non linear space by means of the matrix L, and the 
resulting integration is not univocally defined, simply because the solution S 1-1 = LL' is not unique, since S* -1 
remains unchanged under the substitution L — > LU, where U is an arbitrary unitary matrix. In order to remove this 
arbitrariness, according to Risken[37|. we can work out the integral of the equation of motion in a small time interval 
of length A, by requiring also that the corresponding Markov process: 



a(t n+1 y = a{t n y-A 



S-\t n )d 8 (v-^lnDetS^ (<„ 



k 



(yivi) = Aj = ^s-j(t n ) (22) 



has the correct equilibrium distribution for A — > 0: 



W e q(a) oc VDetS , exp(-^y(a)) (23) 

In fact it is possible to show that, only with the above definition of the drift term, the associated and univocally 
defined Fokker-Planck equation for the probability distribution W(a,t) reads for A — > 0: 



j I i 



(24) 



S~ l ds (v(a)- -LlnDets) 



+ W{a,t) 

which has the equilibrium distribution W eq (a) satisfying: 

]T ^SrldiW eq (a) + W eq (a) ^ .S, fv - ± In Dets\ = (25) 

Indeed, by multiplying both sides of the equations by Skj and summing over j, we obtain the standard equation for 
the equilibrium distribution y/\S\ exp(— (3V). 
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IV. COVARIANT LANGEVIN DYNAMICS FOR IONS AND ELECTRONS 



We want to implement the above formalism in an ab-initio molecular dynamics (MD) at finite temperature dealing 
with electrons and ions within the same formalism, similarly to what was done in the pioneer work by R. Car and 
M. Parrinello[38|]. In the following we will show how the ionic motion can be quite naturally included in the above 
scheme. In fact what we obtained before does not hold only for the electronic parameters, but for a generic set of 
parameters which appear in a variational wavefunction. The ionic positions R can thus be thought as complementary 
parameters. The inclusion of this kind of parameters in the above formalism is straightforward: if M is the number 
of atoms, then S becomes a (p + 3M) x (p + 3M) block-diagonal matrix. The mixed elements Ss a \ mi are always 
zero since wave functions characterized by different sets of atomic positions are orthogonal. Moreover, since the ionic 
positions R belong to the real space, the corresponding metric is the Cartesian one, and is defined by a diagonal 
matrix S(Ri, R r ) — Sn&i,t among all the ion components. We can esplicitly write down the complete set of equations 
for both the atomic and electronic parameters. For the ionic positions we use 

F l (t n , {a(t n }) + x l n 

(X l nX r „) = ^ kr (26) 

with I, r = 1, • • • , 3M and F l being the force acting on the l-ih ionic cartesian coordinate, while for the electronic 
variables Eq. (|22l) holds with i, j = 1, • • • ,p and where — dgV is the force acting on the parameters a, i.e, the gradient 
of the total electronic energy V evaluated at fixed R with respect to these parameters. 

Notice also that the time discretization corresponding to the ionic dynamics is defined by the arbitrary constant Sn 
appearing in the extended metric tensor defined before, namely A^v = AS 1 ^ 1 . It is clear therefore that the relative 
speed between electron and ion dynamics can be tuned to optimize efficiency, exactly as in Car-Parrinello ab-initio 
molecular dynamics. We emphasize here that in the limit A,Ajv — > consistent results are obtained because the 
equilibrium distribution (|23|) remains unaffected by the choice of Sn- 



V. RESULTS AND DISCUSSION 



Once we set up the discretized equations (|22I26|) we can test the above formalism in a simple but realistic case. 
We are going to study the H2 molecule, looking at the temperature behavior of the total energy E and the bond 
distance r between the two hydrogen atoms. We start with this simple system because the above quantities can be 
easily computed, providing therefore useful benchmarks for our technique. 

According to ApplC]the distribution sampled by means of this covariant Langevin dynamics (CLD) represents an 
improvement of the BO only above a temperature T*. At T — our approximate free energy Fq coincides with the 
BO one Fbo, but as soon as T > the Fbo becomes better for T <T* . 

If the temperature is much lower than the electronic gap the BO approximation should be essentially exact and 
can be easily obtained from the potential energy surface (PES) v(r) of the H2 molecule. 

In the following we are going to show that, in this simple system, we cannot distinguish the correct BO low 
temperature behavior and the one implied by our approximate technique, clearly indicating that T* should be almost 
negligible for this system. 

To proceed further we need now to specify what type of correlated variational wavefunction (jj) we adopt in all 
the following calculations, and its dependence on the two electronic positions r\ and r\. In the singlet state the 
orbital function / (71,7=2) is symmetric and positive and is parametrized here as a product of two factors f{fi,?2) = 
/o(ri,r2) x exp( J(fi, ^2)), where fo is taken fixed and allows to satisfy the electron-electron and electron-ion cusp 
conditions, whereas 

J = X>,^; (fiMfa) (27) 

is cusp free and is expanded systematically in a basis of atomic orbitals centered on each atom containing up to 3s 
and lp gaussian functions and a constant one 4>q — 1. This amounts to p — 65 independent variational parameters 
for the symmetric matrix Ajj. The exponent of the gaussians are kept fixed during our simulations. Despite this 
limitation in the choice of the basis this is acceptable for the H2 molecule in a physically relevant range of distances 
between the atoms, as it is shown in Fig.([T]). 

The chosen variational ansatz is particularly useful for evaluating the complicated terms in (|22[) . i.e. the drift- 
diffusion ones which depend linearly from the temperature and require the knowledge of the derivative of the matrix 
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S. This is indeed simpler for the parameters A,,j which appear in a linear fashion in the exponential factor J in 
Eg. ([27)1 . The first step is thus to construct the PES of the molecule ( FigfT]). In this way we not only acquire the key 
information for the numerically exact evaluation of the BO observables, but we also check that our choice of the free 
variational parameters in the wave function allows us to recover the well known PES for this molecule [39l l40j . 




j i i i i i i i i i i i i 

1 1.25 1.5 1.75 2 2.25 2.5 

r [a.u] 



Figure 1: Black line: Total energy E as a function of the bond length r obtained by minimizing the energy of our variational 
wavefunction for fixed r; in doing this we act only on those parameters {a} which are kept free to evolve in the dynamics J22J. 
Red points: Energy with error bars of configurations sampled in a dynamics (|22I26[) with T = 0.01 Ha. The PES is correctly 
followed during the simulation. In the inset a little region around the minimum at r — 1.40 a.u. is enlarged. 



Canonical averages of an observable 0(r) can be obtained by computing numerically the one dimensional (condi- 
tionally convergent) integrals 

q _ fdr r 2 Q(r)exp{-f3v(r)) 
J dr r 2 exp(— /3v(r)) 

On the other hand we can compute O as a time average on the Langevin dynamics (|22I26[) for sufficient low T. The 
extrapolation A — > involving the discretized time steps is performed in the order An —> 0, A — > 0. It is observed 
(see Figj2]) that the An dependence of the time averages of the quantities is linear for fixed A, a property useful in 
the extrapolation. 

Finally we show our results for the total energy and the bond distance at varius temperatures in the range between 
0.001 -j- 0.01 Ha, i.e. from room temperature to ~ 3000 K. The forces acting on the parameters and on the ions, as 
well as the matrix S are evaluated by a short QMC run at each iteration of the dynamics. In Fig.© and in Fig.(HJ) 
we show the outcome of our covariant Langevin dynamics simulations. 

We see that our Langevin dynamics gives result in very good agreement with the expected BO values. We stress once 
again that this dynamics does not require an electronic minimization at each ionic move, realizing an impressive gain 
from the point of view of the computational cost. On the other hand, this kind of dynamics should behave differently 
with respect to the standard BO-MD one when the temperature is raised and for T > T* should be more realistic, 
because corresponding to a more accurate upper bound of the exact free energy F. In figures (|3I4I) we limit the study 
of the average energy and bond length in a range of temperatures smaller than 3000 K because, above this value, 
first dissociation events start t o ap pear during the simulations. This temperature is in good qualitative agreement 
with low pressures experiments [41|. Roughly speaking the dissociation probability depends on the ratio between the 
thermal energy T and the depth of the free energy well AU through the Boltzmann weight [i^] exp(— AU /T) within the 
assumption that excited electronic eigenstates are well-separated in energy from the ground state. There are instead 
examples [22]| in which BO approximation breaks down, particularly near the transition state of a chemical reaction. 
In fact, as the reaction coordinate r increases, the energy gap between the ground state and the first (antibonding) 
excited state becomes smaller |3^, for example when r > 4 a.u. this quantity becomes smaller than 8000if . Therefore 
large fluctuations in the bond length, certainly occurring at large temperatures, are in principle not well described 
under a BO scheme. Since by definition, an atomic dissociation requires to sample correctly events with large r, 
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Figure 2: Time averages of the total energy E at T — 0.003 Ha as a function of An for 4 values of A. All the series converge 
roughly to the same value which is also the expected one (horizontal dashed line) obtained with eq. (|28[) . simplifying the second 
extrapolation A — > 0. 
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Figure 3: Total energy E as a function of temperature. The range of temperature is well below the electronic gap ~ 0.17 Ha 
(see Fig[T]) so the expected exact value is the BO one evaluated by eq. (1281) (black line). Red points are obtained by integrating 
the coupled equations (|22I26I) . Data are in agreement with the predicted values. 



we expect to find differences between the standard BO-MD and the dynamics generated by (|22|26|) . at large enough 
temperatures. In FigJSl we observe that the probability of dissociation is enhanced in our dynamics, which can take 
implicitely into account also the effective repulsion due to the antibonding state. As expected, this is in sharp contrast 
with a DFT-BO dynamics obtained using the Quantum ESPRESSO package [H, [Zj]. In the latter dynamics large 
fluctuations in r do not lead to dissociation, as is partially shown in fig. ([5|). Indeed no escape event occurs within DFT 
BO-MD, even for a long time dynamics. Moreover in order to compensate the well known overbinding error 45] of the 
local density approximation (LDA), we have increased the temperature by a factor proportional to the LDA energy 
barrier (0.2415-ffa), and observed no qualitative changes in the trajectories, always confined around the minimum 
energy value. It is clear therefore that, quite generally, the BO-MD greatly underestimate the evaluation of the 
reaction rate if, for instance, a mean first-passage £ame[42| analysis is performed. 



9 




' 4 0.002 0.004 0.006 0.008 OoT 

T [Ha] 



Figure 4: Bond length r as a function of temperature. Even for this observable the Langevin dynamics (red points) give values 
compatible with the expected ones (black line). 



DFT 




lxlO 4 2xl0 4 3xl0 4 

t [Ha '] 

Figure 5: Bond length r as a function of simulation time at a temperature of T = 8000 K. Coloured points (red, green and 
blue) correspond to simulation performed with the dynamics presented in this work, while the grey solid ones are obtained with 
a DFT- Langevin BOMD. The time step used in the integration of the equations is An = 0.1 Ha -1 and points are plotted 
every 10 iterations. The dashed line indicates the distance r* such that the energy gap between the ground state PES and the 
first excited one becomes smaller than T. All the CLD trajectories show escape events while the DFT one describes a stable 
molecular configuration up to 20 x 10 4 HaT 1 of simulation time (not shown). 



VI. CONCLUSIONS 

In this paper we have introduced a new promising approach to deal with finite temperature simulations of electronic 
systems. The approach is general and, as we have emphasized in the introduction, can be easily extended to sev- 
eral branches of the electronic simulations, from ab-initio finite temperature simulation of realistic systems based on 
Hartree-Fock, DFT or quantum Monte Carlo methods, to finite temperature simulations of strongly correlated Hamil- 
tonians defined on a lattice. In particular this technique allows us to improve systematically the Born-Oppenheimer 
approximation in a temperature range where the quantum effects on atoms are neglegible. In principle also these 
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quantum effects can be dealt in a simple way. To this purpose it is enough to define a quantum ansatz \a) describing 
electrons and ion coordinates quantum mechanically, including in {a} also variational parameters corresponding to 
the atomic wave function $(R), for instance described by gaussians centered around the average atomic positions. In 
that case the same derivation holds as electrons and ions can be dealt in the same footing, the metric matrix S will 
have non trival off diagonal elements between electronic and atomic variational parameters. 

Although our first application is limited to the simple H2 molecule with classical atomic coordinates, this extremely 
simple example already shows that it is possible to catch some qualitatively new behavior, that is not possible to 
describe with the conventional BO approximation. Namely at large enough temperature the molecule can dissociate 
due to non adiabatic effects. 

We plan to extend our method to larger and more complex realistic systems including also quantum effects for 
atoms. Unfortunately, so far we have encountered a difficulty to compute in an efficient way the metric tensor S and 
its derivatives for a generic correlated wave function. For this reason, at present, it looks that the penalty method [34j 
could be a more realistic possibility for extending our technique, because the penalty method does not require the 
evaluation of the derivatives of the metric tensor. Apart for this technical issue there are many open problems that 
can be tackled with this new technique. For instance one would like to know the magnetic transition temperature of a 
piece of material. Without taking into account the electronic entropy this is not possible for most electronic ab-initio 
methods, but, by applying our technique, a reasonable estimate can be easily obtained. In lattice models, an old 
standing problem is for instance the extension of the Gutzwiller variational ansatz to finite temperature calculations. 
Within the variational Monte Carlo it has been established that the Hubbard model for U/t large enough should be 
superconducting with a d-wave order parameter. However it is not possible to predict within the same ansatz the much 
more interesting superconducting temperature and how it depends on the various details of the model, such as doping 
and the value of the Coulomb repulsion U/t. In our formulation what can be done at zero temperature can be readily 
extended to finite temperature and the evaluation of the critical temperature should be straightforward, likewise a 
standard (but much more accurate because including electron correlation) mean field theory at finite temperature. 
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Appendix A: Proof of the integral formula of Eq. (|ll|) 



In this appendix we use known results of differential geometry in Riemann spaces 46] with tensor metric S. In order 
to prove Eq. fTTj) . it is enough to consider the complete basis: 



N 



= IR(n)l > ( A1 ) 



n=l 



where li(n) is an arbitrary choice of N different integers among the L possibilities, that defines the Hilbert of space 
of N fermions containing D = f£) independent states. It is simple to realize that it is enough to prove that, given 
two arbitrary states \i) and \j), we have: 



Oij = j da p ^/\S\{i\a){a\j) = CS id (A2) 

where C is an overall constant. Indeed, by assuming that the above equation holds, we can insert in Eq. (|ll|) the 
completeness I^2t in both the bra and the ket of numerator in Eq. (|lip and obtain: 



/ da p y/\S\(a\ex-p{-PH R )\a) = 
52(j\exp(-pH R )\i)fdaP^}S\(a\j)(i\a) 



CTrexp(-/3ff R ) 

(A3) 



which easily proves Eq. (fTTj) . 
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In order to establish Eq. (|A2|) we can consider the group of transformations a a' that leaves unchanged the metric 
tensor S defined by: 

U\a) = \at) (A4) 

where U is a unitary matrix that maps any variational ansatz a to a new variational ansatz a' of the form defined in 
Eq.Q. To this purpose it is enough to consider the unitary tranformations defined by: 

UiJjj\ - (1 - 2^ m )4, Ui = exp^Trcfd) (A5) 
U P c\ul = c t (J) (A6) 

where p(l) is an arbitrary permutation of the L indices. All the above transformation are real and unitary and 
therefore conserve the distance between two arbitrary vectors, implying that the metric ds 2 remains unchanged under 
all these transformations, when applied to any arbitrary state of the ansatz: 

ds 2 = S lt:j {a)da l da 3 = S hJ (a')da' l da' j (A7) 

In differential geometry these transformations are called isometries, and represent the basis for the classification of 
symmetric Riemann spaces. In this context they are important to prove the main statement of this appendix. Indeed 
we can consider any isometry as a change of variable in the integral and obtain that (since the integration variables 
are dummy variables we can use a in place of a'): 



O id = / da?y/\S\(i\Ul\a)(a\U\j) (A8) 



Now since the set of states is complete the matrix elements Oi j define univocally an operator in the given D— 
dimensional Hilbert space. Therefore by applying the relation (IA8|) for all isometries Ui for I — 1, we obtain 

that this operator O commutes with all fermion occupation number ni and therefore has to be diagonal, namely 
Oij — CiSij. On the other hand we can apply Eq. (|A8[) for an arbitrary unitary permutation Up, that is able to 
connect any state i of the Hilbert space to any other one \ j), namely Up\i) = \j). Thus it easily follows that: 



M = / da?y/\S\(i\U* P \a)(a\U P \i) (A9) 
daPy/\S\{j\a){a\j)=Ojj (A10) 
implying that On = Ci does not depend on i, and this concludes the proof of this appendix. 



Appendix B: Proof of the upper bound for normalized states 



The expectation value of an operator O over a normalized state a is equivalent to average (ipi\0\ipi) over the 
distribution = (ipi\a) 2 over the eigenstates ipi of the operator O. In fact it immediately follows that < pi < 1 and 
that y\. pj = 1. Since for any distribution pi and convex function /, it is well known that, from Jensen's inequality, 
we have: 

(f(H)) > f({H)) (Bl) 

where the symbol (O) means averaging over the distribution pi of the operator O, namely (O) — '^2 i pi{tpi\0\tpi) . Since 
the operator H is Hcrmitian, H and f(H) are diagonalizcd by the same eigenvectors, and therefore the distribution pi 
is the same for both operators and relation (|B1[) simply follows from the convexity of /. Then by using the convexity 
of the function f(x) = exp(— x/T), by applying the above consideration to the operator O = f(H), we obtain: 

(a\ exp(-H/T)\a) > exp((a| - H/T\a)) (B2) 

which concludes the proof of this appendix. 
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Appendix C: Approximate partition function Zq versus exact and Born-Oppenheimer partition functions 

In this appendix we want to investigate the nature of the approximation of the partition function Zq defined in 
(|12[) and used in this work. In order to do this we will compare the approximate partition function Zq with the exact 
Z and the approximate Born-Oppenheimer Z^o-, m the general case when we use p < D variational parameters in 
the normalized wave function ansatz \a). To simplify the notations we avoid to use the dependence on the atomic 
positions R. We assume that the ground state energy E is non degenerate and all the eigenvalues \Ei\ < B, namely 
the spectrum is bounded and B, as well as the maximum gap A = Max^^ — Eq, grows at most linearly with the 
number N of electrons. These assumptions are commonly satisfied by physical Hamiltonians of interacting fermions. 

Within these assumptions, we will see that Zq(T) is an approximation for Z(T) better than Zbo(T) as long as the 
temperature T is larger than a crossover temperature T* < T where T remains finite for N — > oo. 

As mentioned, we assume to know a complete orthonormal set {|*)} i=0 D _i of cigenstatcs of the hamiltonian 
H that operates in a /^-dimensional Hilbert space. This implies that at a given temperature T the exact partition 
function is: 

D-l 



Z(T) = e~ m/T (CI) 



i=0 

whereas the BO partition function is: 

Z BO {T) =exp(-£ y /T) (C2) 

where Ey — Min Q (a|-ff |a) and the approximate partition function Zq is given in Eq. (|12[) . We remind that we have 
already proven, using the convexity of the exponential function, that the relation: 

Z(T) > Z Q (T) (C3) 

holds for every T, and obviously Z(T) > Z BO (T). 

In order to identify the more accurate approximate partition function, namely the one with the larger bound for 
Z(T) we consider the ratio between the Zq and Zbo'- 

^lf> < c4 > 

Since Zq(T) is essentially a classical partition function over p variables, the equipartition theorem immediately implies 
that: 

Z Q {T) ex Z BO (T)TP/ 2 (C5) 

Thus the BO approximation is better at low enough temperature, and, our low temperature free energy Fq — 
Ey — p/2TlnT is expected to be a very bad approximation of the quantum free energy especially when p is very 
large, just because classical and quantum free energy differ substantially at very low temperatures. 

The above consideration could lead to the disappointing conclusion that Cq(^) > 1) namely Fq(T) < Ey, only for 
very high temperatures. 

However we can easily find a lower bound for Cq(^) by using that the spectrum is bounded, as assumed at the 
beginning of this appendix: 



Ca( T) = D !**m<**r*^) > Cexp( _ A/r) (C6) 



When the above bound is larger than one, Cq(T) is certainly larger than one, implying Fq < Fbo- This occurs 
for T > T , where T is easily determined by T = A/lnD. Hence in the thermodynamic limit there exists a finite 
crossover temperature T*, as A/lnD remains finite for N — > oo, according to our assumptions. 



Appendix D: Slater determinants and symmetric Riemann spaces 



We consider the space A4 of normalized Slater determinants in a finite dimensional Hilbert space H where fermions 
can occupy L different one particle states, denoted by conventional creation operators c\. A Slater determinant with 
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N electrons can be formally written in second quantization notations by means of N x L real numbers ip)'- 



N L 



(Dl) 



i=ij=i 



However all the variables of the matrix ip are highly redundant because, as well known, the Slater determinant after 
the linear transformations ip ~~ ► Lip is multiplied by a constant \ip) — ¥ \L\\ip), whereL is an arbitrary N x N matrix 
and \L\ its determinant. It is clear that, in order to define a Slater determinant with unit norm we can consider 
one constraint (ip\ip) — \ipip >\ — 1 over the NL variables defining the N x L matrix ip, amounting therefore to 
NL — 1 independent real variables. By the above discussion, the wavefunction \tp) is left invariant for all matrix 
transformation ip — > Lip with \L\ = 1, defining N 2 — 1 independent variables for L. Thus, it follows that \ip) can be 
parametrized by (NL — 1) — (N 2 — 1) = N(L — N) independent real variables. In a more rigourous mathematical 
formalism, by neglecting an immaterial overall sign ±1 in the definition of ip, the space M. represents the coset space 
0(L,L- N), where 0(L,L- N) is the irreducible symmetric Riemannian space SO{L)/S(0{N) x 0(L - JV)).{46| 
We remind here that O(N) denotes the group of generic orthogonal matrices, whereas SO(N) represents the group 
of orthogonal matrices with determinant one. Similarly O(N) x 0(L — N) represents the group of block diagonal 
matrices with N x N and L — N x L — N blocks, where each block is in turn an orthogonal matrix. Also the symbol 
S(0(N) x 0(L — N)) indicates that the determinant of this block diagonal matrix (the products of the determinant 
of each block, equal to ±1 as for any orthogonal matrix) has to be 1. 

This space M. is compact (all the N(L — N) independent variables represent essentially angles of unit vectors in L 
dimensional space) and there exist a unique (up to a constant) measure d\i such that dUfi = d\x for all U G SO(L) 
where SO(L) is the group of L x L orthogonal matrices with unit determinant [46j . namely \U\ = 1. An orthogonal 
matrix U, acts on \ip) in an obvious way, namely ip — > ipU in Eq. (|Dip . The space M. can be therefore represented 
by an irreducible symmetric Riemannian space. Using a matrix U G SO(L) we have essentially L orthonormal 
directions (e.g. the raws of the matrix), and the hrst N spans all possible Slater determinants in the space M.. For 
the previous discussion this Slater determinant will be left unchanged (up to a sign) if we multiply the matrix U 
for an arbitrary element of the S(0(N) x 0(L — N)) unitary group, and therefore M. is equivalent to the space 
SO(L)/S(0(N) x 0(L — N)), 

As a further proof that M. is equivalent to SO(L)/S(0(N) x 0(L — N)) it is also easy to verify that the dimension 
of this space space is exactly N(L — N). The dimension of an orthogonal matrix of dimension D is D(D — l)/2, and 
therefore the dimension of the coset space SO(L)/S{0(N) x 0(L - N)) is L(L - l)/2 — (L — N)(L — N — l)/2 — 
N(N — l)/2 = N(L — N)M. 

In order to represent the irreducible space SO(L)/S(0(N) x 0(L — N)) for L » N, with N(L — N) variables, a 
possible choice is to define an unconstrained N x (L — N) matrix V and the corresponding unitary L x L matrix Q: 



with the constraint that the positive definite matrix VV^ has all eigenvalues bounded by one, namely VV^ < 1. Thus 
we explicitly see that the space is compact. As mentioned before we can identify a wavefunction ip G M with the first 
N raws of this unitary matrix Q, up to a sign, so that the orbitals of the determinant are: 



A measure dip of the coset (reducible) Riemann space SO(L) /S(0(N) x 0(L — N)) is said to be an invariant measure 
when it remains invariant under all unitary transformations U G U(L), namely dipU — dip. An invariant measure 
represented by the matrix V is given by: 



where C is an appropriate normalization constant, and fi(V) is the invariant measure in SU(L)/S(U(N) x U(L — 
AO). [461] Although explicit formulas are known for the invariant measure, they look a bit complicated to be implemented 
in practice. We are confident that a very convenient expression of the invariant measure is possible in terms of the 
eigenvalues of W^, which should amount to only TV 3 operations. This would lead immediately to a computationally 
affordable extension of our method to DFT or mean-field type of ansatz. 




(D2) 



ipi,k = Qi,k for 1= 1,2,- •• ,N. 



(D3) 



dip = Cd/i(V) 



(D4) 
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